ChIP AML PiPeline v2¶
In [2]:
import os
from IPython.display import IFrame
import pandas as pd
import sys
sys.path.insert(0, '../..')
from JKBio.epigenetics import ChIP_helper as chiphelper
from JKBio import Helper as helper
import igv
import numpy as np
import pyBigWig
import itertools
import seaborn as sns
from scipy import stats
from matplotlib import cm
from scipy.cluster.hierarchy import linkage, leaves_list
from matplotlib import pyplot as plt
from sklearn.cluster import AgglomerativeClustering
from sklearn.mixture import GaussianMixture
from bokeh.plotting import *
from sklearn.manifold import MDS, TSNE
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale
from pybedtools import BedTool
output_notebook()
%load_ext autoreload
%autoreload 2
adding the data bucket to path¶
In [ ]:
! gcsfuse --only-dir Chip_AML jkobject data/seqs
processing using Nextflow¶
In [ ]:
singleend, pairedend = chiphelper.extractPairedSingleEndFrom('data/seqs')
Pipeline¶
- Raw read QC (FastQC)
- Adapter trimming (Trim Galore!)
- Alignment (BWA)
- Mark duplicates (picard)
- Merge alignments from multiple libraries of the same sample (picard)
- Re-mark duplicates (picard)
- Filtering to remove: blacklisted regions, duplicates, primary alignments,unmapped,multiple locations, containing > 4 mismatches, insert size > 2kb, map to different chromosomes
- Alignment-level QC and estimation of library complexity (picard, Preseq)
- Create normalised bigWig files scaled to 1 million mapped reads (BEDTools, bedGraphToBigWig)
- Generate gene-body meta-profile from bigWig files (deepTools)
- Calculate genome-wide IP enrichment relative to control (deepTools)
- Calculate strand cross-correlation peak and ChIP-seq quality measures including NSC and RSC (phantompeakqualtools)
- Call broad/narrow peaks (MACS2)
- Annotate peaks relative to gene features (HOMER)
- Create consensus peakset across all samples and create tabular file to aid in the filtering of the data (BEDTools)
- Count reads in consensus peaks (featureCounts)
In [ ]:
! nextflow cloud create 'JKcluster' -c 4
In [ ]:
! nextflow cloud create jkcluster -c "nextflow/nextflow.config" 40 && \
nextflow nf-core/chipseq -c "nextflow/nextflow.config" \
--singleEnd \
--seq_center 'DFCI' \
--email 'jkobject@gmail.com' \
--bucket-dir 'gs://jkobject/Chip_AML/nextflow/CHIPprocess_2/' \
--keyfile '~/jkobject-b6f1adaffcb8.json' \
--projectname 'jkobject' \
--zone 'us-east1-b' \
--skipDiffAnalysis \
--narrowPeak \
--design "nextflow/design.csv" \
--genome 'GRCh38' \
--profile gcp \
--resume \
--skipPreseq \
--max_cpus 8 && \
nextflow cloud shutdown jkclustert
Gathering data¶
In [ ]:
!gsutil -m cp -r gs://amlproject/Chip/results/bwa/mergedLibrary/macs/narrowPeak/ ../data
In [ ]:
!gsutil -m cp -r gs://amlproject/Chip/results/bwa/mergedLibrary/bigwig/ ../../data
In [ ]:
!cp ../../data/narrowPeak/*MV411*.narrowPeak ../../data/MV4narrow
In [110]:
! mkdir ../../data/BroadPeaks/MV411 && mv ../../data/BroadPeaks/MV411_* ../../data/BroadPeaks/MV411/
In [26]:
bindings = chiphelper.loadPeaks('../../data/MV4narrow/', isMacs=False,skiprows=0)
In [27]:
broadbindings = chiphelper.loadPeaks('../../data/BroadPeaks/MV411/', isMacs=False,skiprows=0)
In [28]:
SEgenes = pd.read_csv('../data/SEgenes.csv')
CTF = pd.read_csv('../data/CTF.csv', header=None)[0].tolist()
In [29]:
CTF.extend(['GATA2','IKZF1','LYL1' ,'PU1','SMC1'])
CTF
Out[29]:
In [31]:
CTF = list(set(CTF))
In [32]:
peaks = !ls ../../data/MV4narrow/*.narrowPeak
broadpeaks = ! ls ../../data/BroadPeaks/MV411/*.broadPeak
peaks = set([i.split('/')[-1].split('.')[0] for i in broadpeaks]) | set([i.split('/')[-1].split('.')[0] for i in peaks])
looking at the data and renaming¶
In [33]:
peaks
Out[33]:
In [34]:
set(bindings['name'])
Out[34]:
In [35]:
bindings
Out[35]:
In [36]:
broadbindings
Out[36]:
In [38]:
bindings = bindings[~bindings.name.isin(set(broadbindings.name))]
In [40]:
bindings = bindings.append(broadbindings)
In [41]:
len(bindings)
Out[41]:
In [42]:
bindings['replicate']= [i.split('-')[-1][-1] for i in bindings['name']]
bindings['tf'] = [i.split('-')[2] for i in bindings['name']]
In [43]:
bindings['peak_number'] = ['_'.join([i.split('_')[2],i.split('_')[5]]) for i in bindings['peak_number']]
In [44]:
bindings
Out[44]:
In [45]:
bindings.to_csv('../temp/bindings.bed',sep='\t',index=False)
In [3]:
bindings= pd.read_csv('../temp/bindings.bed',sep='\t',header=None, index_col=None,
names=["-log10pvalue","-log10qvalue", "chrom", "end", "foldchange", "name", "peak_number", "relative_summit_pos", "start", "replicate","tf"])
In [276]:
bindings
Out[276]:
In [20]:
from gsheets import Sheets
sheets = Sheets.from_files('~/.client_secret.json', '~/.storage.json')
url="https://docs.google.com/spreadsheets/d/1yFLjYB1McU530JnLgL0QIMAKIkVl3kl0_LCHje2gk8U"
gsheet = sheets.get(url).sheets[2].to_frame()
gsheet
Out[20]:
In [21]:
bw = ! ls ../../data/bigwig
bw
Out[21]:
In [46]:
len(set(bindings.name))
Out[46]:
In [25]:
len(bw)
Out[25]:
In [48]:
# ONE off
for i in bw[2:]:
a = gsheet[gsheet.id=='mp'+i.split('_')[2]].name.values[0]
i = '../../data/bigwig/'+i
a = '../../data/bigwig/'+a+'.mLb.clN.bigWig'
! mv $i $a
print(a)
In [22]:
set(bindings.name)
Out[22]:
In [ ]:
replicates = chiphelper.findReplicates(folder='data/seqs/results/bwa/', sep='_', namings='_R([0-9])',namepos=0)
Visual inspection of the features and and look at QCs¶
process:¶
look at all t with a very low frip score as noted by encode.
look at all peaks tracks together and see for location of intense co binding.
- if we can discern peaks and if, for some reasons, some good peaks are not called by macs.
- if looks good and we can see a lot of peaks.
- if a lot of noise but seems consistent with replicates.
- if just seems to have very few peaks.
Validate still but flag as potentially bad.
Else remove.
results:¶
In [4]:
bad=[
"mp168",
"mp129",
"mp128",
"mp773",
"mp774",
"mp575",
"mp614",
"mp714",
"mp433",
"mp156",
"mp650",
"mp604",
"mp27",
"mp627",
"mp117",
"mp771",
"mp118",
"mp431",
"mp430",
"mp324",
"mp565",
"mp569",
"mp125",
"mp627",
"mp568",
"mp427",
"mp124",
"mp716",
"mp581",
"mp589",
"mp321",
"mp601",
"mp745",
"mp772",
"mp770",
"mp590",
"mp623",
"mp718"]
merging duplicates¶
In [6]:
%matplotlib inline
mergedpeak, tomergebam, remove, ratiosofunique = chiphelper.mergeReplicatePeaks(bindings,'../../data/bigwig/',markedasbad=bad, window=150, mincov=4, doPlot=True, minKL=10, cov={}, use='poisson', MINOVERLAP=0.25,lookeverywhere=True, only='',saveloc='')
In [ ]:
tomergebam
In [ ]:
mergedpeak = mergedpeak[mergedpeak.columns[[2,9,3,5,6,4,0,1,7,10]]]
In [ ]:
mergedpeak.to_csv('../temp/merged_replicates.csv')
In [7]:
mergedpeak = pd.read_csv('../temp/merged_replicates.csv', index_col=0)
show replicates overlap¶
sorting and removing samples¶
In [7]:
bigwigs=os.listdir('../../data/bigwig/')
for val in bigwigs:
for v in remove + toremove + ['scale','POLII','IGG','CTCF','INPUT']:
if v in val:
bigwigs.remove(val)
break
bigwigs = ['data/bigwig/'+ i for i in bigwigs]
In [62]:
set(mergedpeak.tf)
Out[62]:
In [63]:
mergedpeak.foldchange.min()
Out[63]:
In [ ]:
mergedpeak['name']=mergedpeak.tf
In [ ]:
## Removing bad ChIP protein
mergedpeak = mergedpeak[~mergedpeak['name'].isin(['CDK13','GSE1'])]
Create a consensus set¶
In [ ]:
merged = chiphelper.simpleMergedPeaks(mergedpeak[~mergedpeak.tf.isin(['MED1','SMC1','CTCF','POLII','IRF2BP2_FLAG','IRF2BP2', 'H3K27ac', 'H3K27me3', 'H3K4me3', 'H3K79me2',])], window=150)
In [49]:
len(merged)
Out[49]:
In [8]:
len(mergedpeak)
Out[8]:
In [53]:
merged
Out[53]:
In [ ]:
merged.to_csv('../temp/merged.bed', sep='\t',index=None)
In [11]:
merged = pd.read_csv('../temp/merged.bed', sep='\t')
In [13]:
fig = sns.pairplot(merged[merged.columns[8:14]], corner=True, diag_kind="kde", kind="reg", plot_kws ={"scatter_kws":{"alpha":.05}})
def col_nan_scatter(x,y, **kwargs):
df = pd.DataFrame({'x':x[:],'y':y[:]})
df = df[df.sum(0)!=0]
x = df['x']
y = df['y']
plt.gca()
plt.scatter(x,y)
def col_nan_kde_histo(x, **kwargs):
df = pd.DataFrame({'x':x[:]})
df = df[df['x']!=0]
x = df['x']
plt.gca()
sns.kdeplot(x)
fig = fig.map_upper(col_nan_scatter)
fig = fig.map_upper(col_nan_kde_histo)
fig.savefig('../results/cobinding/pairplot_experiments.pdf')
plt.show()
counts,val = np.unique(merged[merged.columns[8:]].astype(bool).sum(1).values, return_counts=True)
fig = sns.barplot(data=pd.DataFrame(val, index=counts,columns=['counts']).T).set_yscale("log")
fig.savefig('../results/cobinding/pairplot_experiments.pdf')
plt.show()
i = merged[merged.columns[8:]].astype(bool).sum(1)
print(i.max(),i.mean(),i.min())
In [14]:
counts,val = np.unique(merged[merged.columns[8:]].astype(bool).sum(1).values, return_counts=True)
fig = sns.barplot(data=pd.DataFrame(val, index=counts,columns=['counts']).T).set_yscale("log")
TODO: what would be the distribution in case they are completly random?
Correlation over consensus set¶
In [38]:
sns.clustermap(np.corrcoef(stats.zscore(merged[merged.columns[8:]].values.T, axis=1)), figsize=(20, 20), xticklabels=merged.columns[8:], yticklabels=merged.columns[8:]).ax_col_dendrogram.set_visible(False)
In [445]:
merged[merged.columns[:8]].to_csv('../temp/conscensus.bed',sep='\t',index=None, columns=None)
In [ ]:
cor = np.log2(1+cor)
a = cor.var(0).argsort()
subcor = cor[:,a[-TOTVAR:]]
In [ ]:
subcor.shape
In [ ]:
subcor = PCA(PCS).fit_transform(subcor)
In [ ]:
%matplotlib inline
cmaps = ['Greys', 'Purples', 'Blues', 'Greens', 'Oranges', 'Reds',
'YlOrBr', 'YlOrRd', 'OrRd', 'PuRd', 'RdPu', 'BuPu',
'GnBu', 'PuBu', 'YlGnBu', 'PuBuGn', 'BuGn', 'YlGn']
import seaborn as sns
from matplotlib import pyplot
fig, ax = pyplot.subplots(figsize=(20,15))
sns.heatmap(ax=ax, data=subcor,xticklabels=False, cmap=cmaps[-4],
yticklabels=[i.split('.')[-4].split('/')[-1] for i in bigwigs],
cbar_kws={"orientation": "horizontal"},
vmax=0.7).set_title('PCs of each binding profile over conscensus peakset')
#fig.savefig("temp/co_occupancy_matrix.png")
In [ ]:
fig.savefig('PCs_binding_peakset.png')
clustering¶
I have tried gaussian mixtures and Agglomerative clustering algorithm. Only the second can create a hierarchical clustering.
It seems that gaussian mixture makes more sense given the data we have, for now, is more "homogeneous".
I am still not so happy with the clustering. It can be because of the how much importance, outlier values and the high number of noisy values from locations with no peaks.
We can use similar methods to RNAseq to improve this (clamping values, log transform, first round of PCA..)
In [ ]:
labels = GaussianMixture(n_components=2, covariance_type='diag').fit_predict(subcor)
In [ ]:
names = np.array([i.split('.')[-4].split('/')[-1] for i in bigwigs])
sort = labels.argsort()
p = helper.plotCorrelationMatrix(data=cor[sort],
names=names[sort],
colors=labels[sort],
title="correlation between TF occupancy",
interactive=True)
In [ ]:
show(p)
In [ ]:
p = helper.scatter(TSNE(2,5).fit_transform(subcor),labels=names, colors=labels)
In [ ]:
show(p)
In [ ]:
sns.clustermap(subcor)
annotatePeaks¶
In [5]:
additional = {}
additional['activation'] = chiphelper.simpleMergedPeaks(mergedpeak[mergedpeak.tf.isin(["H3K27ac",'H3K79me2','H3K36me3','H3K4me3'])], window=10, mergedFold="max")
additional['repression'] = mergedpeak[mergedpeak.tf=='H3K27me3']
additional['IRF2BP2'] = mergedpeak[mergedpeak.tf=='IRF2BP2_FLAG']
additional['MED1'] = mergedpeak[mergedpeak.tf=='MED1']
additional['SMC1'] = mergedpeak[mergedpeak.tf=='SMC1']
additional['CTCF'] = mergedpeak[mergedpeak.tf=='CTCF']
additional['POLII'] = mergedpeak[mergedpeak.tf=='POLII']
In [ ]:
#add enhancers at TSS in the matrix (promoters)
# add super enhancers
# add other enhancers
In [ ]:
# liftover enhancer assignements
# merge enhancer assigned to cobinding
In [ ]:
chiphelper.MakeSuperEnhancers()
In [6]:
for key, val in additional.items():
merged[key] = chiphelper.putInConscensus(merged[merged.columns[:8]],val)
merged = merged.replace(np.nan,0)
In [7]:
merged = merged.replace(np.nan,0)
In [8]:
merged.to_csv('../temp/merged.bed', sep='\t',index=None)
In [3]:
merged = pd.read_csv('../temp/merged.bed', sep='\t')
Co Binding Matrix¶
Look at AUC for all ChIPs over all peaks of all ChIPs
In [178]:
rand = np.random.choice(merged.index,5000)
viridis = cm.get_cmap('viridis', 256)
data = merged[merged.columns[-7:]]
for val in data.columns:
data[val] =stats.zscore(np.log2(1+data[val]))
data = (((data -data.min(0))/ (data.max(0)))*256).astype(int)
data = data.loc[rand]
for val in data.columns:
a = [viridis(v) for v in data[val]]
data[val] = a
In [185]:
fig = sns.clustermap(np.log2(1.01+merged[merged.columns[8:-7]].loc[rand].T),col_cluster=False, z_score=0, vmin=0,vmax=3, col_colors = data, figsize=(30,20),xticklabels=False)
fig.ax_col_dendrogram.set_visible(False)
fig.savefig('../temp/clustermap_cobinding_scaled_marks.pdf')
plt.show()
In [198]:
fig = sns.clustermap(np.log2(1.01+merged[merged.columns[8:-7]].loc[rand].T), vmin=0,vmax=3,figsize=(30,20),z_score=0,col_colors=data, xticklabels=False)
fig.ax_col_dendrogram.set_visible(False)
fig.savefig('../temp/clustermap_cobinding_scaled_marks.pdf')
plt.show()
Looking at peak overlap¶
How many of peak in A (column) overlaps with peak in B (rows)
in other words:
what is the percentage of B's peaks that are overlaped by A's peaks
In [90]:
correlation, overlap = chiphelper.computeOverlapCorrelation(merged, norm=True)
In [274]:
data = pd.DataFrame(data=overlap,index=merged.columns[8:], columns=merged.columns[8:])
link = linkage(data.iloc[:-7]) # D being the measurement
col = data[-7:]
col = col[[co for co in col.columns if co not in col.index.tolist()]]
for val in col.columns:
a = [viridis(v) for v in col[val]]
col[val] = a
fig = sns.clustermap(data.iloc[:-7][data.columns[np.concatenate((leaves_list(link),[26,27,28,29,30,31,32]))]], row_linkage=link, col_colors=col.T, col_cluster=False)
fig.ax_col_dendrogram.set_visible(False)
fig.savefig('../results/cobinding/overlap.pdf')
plt.show()
data on the experiments¶
In [193]:
pd.concat([merged[merged.columns[8:]].astype(bool).sum(0),
merged[merged.columns[8:]].max(),
merged[merged.columns[8:]].replace(0, np.NaN).mean(),
merged[merged.columns[8:]].replace(0, np.NaN).var()],axis=1).rename(columns={0:'sum',1:'max',2:'mean',3:'std'})
Out[193]:
Correlation only on overlaps¶
on the overlaps given above
In [273]:
data = pd.DataFrame(data=correlation,index=merged.columns[8:], columns=merged.columns[8:])
link = linkage(data.iloc[:-7]) # D being the measurement
col = data.iloc[-7:]
col = col[[co for co in col.columns if co not in col.index.tolist()]]
for val in col.columns:
a = [viridis(v) for v in col[val]]
col[val] = a
fig = sns.clustermap(data.iloc[:-7][data.columns[np.concatenate((leaves_list(link),[26,27,28,29,30,31,32]))]], row_linkage=link, col_colors=col.T, col_cluster=False)
fig.ax_col_dendrogram.set_visible(False)
fig.savefig("../results/cobinding/correlation_onoverlap.pdf")
plt.show()
In [ ]:
chiphelper.filterOnEnrichment(,toremove=[""])
TODO¶
plots during merge¶
- plot a histogram (bincount) of the signal for each peaks:
- peaks in A, B
- unique peaks in A, B
- peaks in merged version (avg)
- scatter plot matrix of signal correlation between 1to1 for all
- print average size in unique vs overlap peaks
- ratios of unique peaks to total peaks for each TF
plots after merge¶
- create a correlation plot on 1to1 overlap only.
- histogram of number of peaks per amount of TF cobound
- create a correlation matrix only on A's peakset
- create a normalized overlap matrix (given the number of peaks and total peakset, what is the expected overlap and how far away are we from it) z-score norm
- create a clustermap over peaks
improving merge¶
- analyse only when peak is supported by two replicates
- remove GSE1 and CDK13 from cobinding matrix
- do quantile normalization over the signals (after setting up the zeroes)
- removing bad chips:
- look, when badchip only does not correlate with anything else..
- do log2+1 transform
peak enrichment¶
- use chromHMM
- use the Chipseq data we have (RNApol2, CTCF...)
- filter weird enrichments (not in H3k27ac? in H4me3..)
- make super enhancer enrichment from superenhancer matrix
- merge the clusters and display enrichment of each (in TF and other marks)
Next step:¶
- do cobinding based on gene targeted from ABC-model and merge the same way as before (mean value); merge enrichments as well
- look at data targetting CRC members
- look for enrichment at the gene level (from gene annotations)
- compute/refine binding patterns
- use MEME
- use DeepBind
- look for difference in binding patterns to the litterature (db of binding patterns)
- look for predictability of RNP expression change based on bound genes
- do we see genes with opposite effects?
- can it be explained by the set of co-bound TFs?
- do dependency/expression prediction from all gene level features
- look for most explanatory regressors
- repeat the process across all samples with H3K27ac+RNAseq data we have.
- call mutations from H3k27ac data
- MEME analysis of likely bound TFs, except if TF is not expressed
- compute enrichments
- from same gene assignements, as found on the general consensus peak set, can we find good dependency/expression prediction
- does location
questions:¶
- how did andrew defined significant peaks?
based on closest expressed gene¶
In [ ]:
chiphelper.AssignToClosestExpressed()
based on the ABC model¶
They tested a new model based on and validated by CRISPRi-FlowFISH which is basically able to find enhancer mapping to genes. They used it to compute their model's Accuracy and found a 70% accuracy compared to less than 50% for closest expressed gene.
Way to integrate our HiC data (need ATAC-seq like data as well, but openly available)
In [ ]:
Helper.scatter(TSNE(2,5).fit_transform(data.T), labels=zones.columns[11:],colors=labels)
Compute the CRC and merge with gene assignement¶
~1500 lines of code. Seems to be a slightly updated version from the code I used the first time which was from their originial CRC paper.
There is not a lot of documentation but it was just updated last week and might continue like that
find set of TFs that explain the most cobindings¶
relate to gene depencies¶
what else?¶
Compare data with other labs (H3K27, HiC..)¶
we need to redo everything for similar normal cell type, getting TFs based on the CRC (find it with CRCmapper or on litterature)